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I .  INTRODUCTION 


In  this  report  we  will  be  concerned  with  studying  the  propagation 
and  interaction  of  solitary  waves  in  a face-centered-cubic  (FCC) 
lattice.  Our  intention  will  be  to  determine  the  extent  to  which  the 
properties  of  these  waves,  well  known  in  many  one-dimensional  models, 
are  affected  by  the  additional  space  dimensions  and  to  assess  the  likely 
physical  implications  of  these  properties.  Most  of  the  work  will  be 
undertaken  using  a computer-molecular-dynamic  model,  although  extensive 
analytic  work  is  also  presented. 


Since  the  discovery  that  solitary-wave  solutions  of  the  Korteweg- 
de  Vries  (KdV)  equation  were  stable  to  mutual  collisions,  these  stable 
pulses  (solitons)  have  received  considerable  attention.  Interest  in  the 
problem  arises  from  the  fact  that,  if  relatively  stable  pulses  exist  in 
real  physical  systems,  they  will  substantially  affect  the  manner  in 
which  energy  is  transported  in  the  system  as  well  as  the  speed  with 

2 

which  it  approaches  thermal  equilibrium.  Early  work  on  solitons  was 
devoted  to  the  study  of  one-dimensional  systems  whose  equations  of 
motion  could  at  least  be  approximated  by  equations  which  possessed 

soliton  solutions.  More  recently,  some  progress3  6 has  been  made  in 
extending  the  results  to  multiple  dimensions,  although  the  work  is 
largely  mathematical  in  nature  and/or  applicable  to  only  rather  idealized 

7 

systems.  In  more  applied  work,  Schneider,  Stoll,  and  Hiwatari  have 
observed  in  a three-dimensional  discrete- lattice  model  that  heat  pulses 

exhibited  solitonlike  properties,  and  Ikezi  has  experimentally  verified 
the  existence  of  solitons  in  plasmas.  These  latter  results  suggest  that 


N.J.  Zabusky  and  M.D.  Kruskal,  "Interaction  of  Solitons  in  a Col- 
lisionless Plasma  and  the  Recurrence  of  Initial  States,"  Phys. 

Rev.  Letters  Jj^,  240  (1965). 

A.C.  Scott,  F.Y.F.  Chu,  and  D.W.  McLaughlin,  "The  Soliton:  A New 
Concept  in  Applied  Science,"  Proc.  IEEE  61_,  1443  (1973). 

J.  Denavit,  N.R.  Pereira,  and  R.N.  Sudan,  "Two-Dimensional  Stabil- 
ity of  Langmuir  Solitons,"  Phys.  Rev.  Letters  33,  1435  (1974). 

K. H.  Spatschek,  P.K.  Shukla,  and  M.Y.  Yu,  "On  the  Two-Dimensional 
Stability  of  Ion-Acoustic  Solitons,"  Phys.  Letters  54A,  419  (1975). 

Y.  Gell,  "Drift  Solitons  and  Their  Two-Dimensional  Stability," 

Phys.  Rev.  A 16,  402  (1977). 

R.G.  Newton,  "Three-Dimensional  Solitons,"  J.  Math.  Phys.  lj>,  1068 
(1977) . 

T.  Schneider,  E.  Stoll,  and  Y.  Hiwatari,  "Solitonlike  Properties 
of  Heat  Pulses,"  Phys.  Rev.  Letters  39,  1382  (1977). 

H.  Ikezi,  "Experiments  on  Solitons  in  Plasmas,"  in  Solitons  in  Action 
edited  by  K.  Lonngren  and  A.  Scott  (Academic,  New  York,  1978), 
p.  153. 


soliton  propagation  may  indeed  be  an  important  real  effect  in  solids 
and  further  studies  are  appropriate  at  this  time. 

Our  motivation  for  investigating  the  particular  problem  of  solitary 
waves  in  the  FCC  lattice  arose  from  our  efforts  to  explain  anomalous 
effects  observed  in  computer  simulations  of  shock  waves  in  discrete 
lattices.  The  first  studies  of  this  problem  were  the  three-dimensional, 
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computer-molecular-dynamic  calculations  reported  by  Tsai  and  his  co- 
workers approximately  ten  years  ago.  Their  work  suggested  that,  though 
the  shock  front  propagates  at  a constant  velocity,  the  shock  profile 
itself  is  not  steady.  In  fact,  the  transition  zone  which  connects  the 
two  equilibrated  regions  ahead  of  and  behind  the  front  appeared  to  in- 
crease in  length  as  the  shock  wave  propagated  farther  into  the  lattice. 

In  addition,  the  energy  associated  with  random  motion  appeared  to  be 
higher  directly  behind  the  shock  front  than  in  the  equilibrated  region 

far  behind  the  front.  Shortly  thereafter  Tasi^  ^ reported  the  results 
of  his  studies  of  shock  propagation  in  a one-dimensional,  quiescent 
lattice.  He  also  observed  a nonsteady  shock  profile  and  explained  its 
occurrence  on  the  basis  of  solitary-wave  propagation. 

The  shock  profile  predicted  by  these  computer  simulations  contra- 
dicted the  shock-wave  structure  predicted  by  the  usual  continuum 
equations  and,  if  true,  could  have  a significant  impact  on  several  areas 
of  ballistics  research.  In  view  of  this  fact,  we  initiated  an  in-house, 
computer-molecular-dynamics  program  to  study  shock  waves  in  crystal 
lattices.  We  initially  restricted  our  investigations  to  studies  of 
shock  structure  in  a one-dimensional  lattice.  Our  calculations  in  the 
quiescent  lattice  essentially  confirmed  Tasi's  earlier  results.  In 
addition,  however,  we  found  that  stable  solitary  waves  were  also  gener- 

13  14 

ated  when  the  lattice  was  initially  at  some  nonzero  temperature. 

9i  D.H.  Tsai,  "An  Atomistic  Theory  of  Shock  Compression  of  a Perfect 

Crystalline  Solid,"  in  Accurate  Characterization  of  the  High-Pressure 
Environment,  edited  by  E.C.  Lloyd,  Natl.  Bur.  Stds.  Spec.  Publ. 

No.  326  (U.S.  GPO,  Washington,  DC,  1971),  p.  105. 

10.  J.  Tasi,  "Perturbation  Solution  for  Growth  of  Nonlinear  Shock  Waves 
in  a Lattice,"  J.  Appl.  Phys.  43^,  4016  (1972).  See  also  Erratum  [J . 
Appl . Phys.  44,  1414  (1973)]. 

11.  J.  Tasi,  "Far-Field  Analysis  of  Nonlinear  Shock  Waves  in  a Lattice," 
J.  Appl.  Phys.  44,  4569  (1973). 

12.  J.  Tasi,  "Perturbation  Solution  for  Shock  Waves  in  a Dissipative 
Lattice,"  J.  Appl.  Phys.  4£,  2245  (1973). 

13.  J.H.  Batteh  and  J.D.  Powell,  "Shock  Propagation  in  the  One-Dimen- 
sional Lattice  at  a Nonzero  Initial  Temperature,"  J.  Appl.  Phys. 

49,  3933  (1978). 

14.  J.H.  Batteh  and  J.D.  Powell,  "Soliton  Propagation  in  a One-Dimen- 
sional Lattice  Under  Shock  Compression,"  in  Solitons  in  Action, 
edited  by  K.  Lonngren  and  A.  Scott  (Academic,  New  York,  1978)  p.  257. 
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In  an  effort  to  determine  whether  the  unexpected  results  observed  in 
three  dimensions  could  be  explained  on  the  basis  of  solitary-wave 
behavior,  we  modified  our  code  to  treat  shock  propagation  in  a face- 
centered-cubic  lattice.  The  results  of  that  investigation  will  be 
presented  in  a subsequent  report.  Here,  we  mention  only  that  the  sim- 
ulations indicated  the  presence,  at  least  for  early  times,  of  solitary 
waves  near  the  shock  front.  It  soon  became  apparent  that,  in  order  to 
interpret  the  results  of  the  shock-wave  calculation,  we  would  first 
have  to  understand  the  behavior  of  solitary  waves  in  three-dimensional 
lattices.  Consequently,  we  initiated  the  studies  reported  here  which 
were  designed  solely  to  elucidate  the  properties  of  solitary  waves  in  a 
three-dimensional  system.  Specific  points  we  include  in  this  report 
include  the  extent  to  which  the  FCC  lattice  can  support  solitary  waves; 
the  effect  of  mutual  collisions  upon  the  wave  profile;  the  stability 
of  the  profiles  to  perturbations  in  the  direction  of  propagation  as 
well  as  perpendicular  to  it;  and  the  extent  to  which  longitudinal  and 
transverse  oscillations  are  coupled  and  how  the  energy  is  exchanged 
between  them. 

We  begin  in  Sec.  II  by  describing  the  model  under  consideration, 
writing  down  the  equations  of  motion  for  each  atom  in  the  lattice,  and 
describing  the  method  for  solving  them.  In  Sec.  Ill  we  present  the 
results  of  the  numerical  studies  which  have  been  undertaken  using  the 
discrete-lattice  model.  In  particular,  we  address  the  generation  of 
solitary  waves  and  discuss  their  stability.  It  is  pointed  out  that  in 
some  cases  it  is  possible  to  generate  coupled  longitudinal  and  transverse 
solitary  waves  which  propagate  in  phase  at  the  same  velocity.  In  Sec. 

IV  we  derive  the  long-wavelength  (continuum)  limit  of  the  equations  of 
motion  for  the  FCC  lattice.  It  is  demonstrated  that  there  exists  a 
steady  solution  to  the  equations  which  predicts  the  coupled  solitary 
waves  found  in  the  previous  section.  The  equations  are  solved  numeric- 
ally and  an  approximate  analytic  solution  is  also  obtained.  Finally, 

Sec.  V contains  a summary  of  the  results  and  the  conclusions  drawn  as 
well  as  some  discussion  of  the  physical  implications  that  the  existence 
of  solitary  waves  might  have  on  energy  transfer  in  crystalline  solids. 


II.  MODEL  AND  EQUATIONS  OF  MOTION 


The  model  whose  properties  we  wish  to  study  consists  of  a pure, 
monatomic,  FCC  lattice  which  is  made  as  long  as  necessary  in  the  z 
direction  and  which  is  periodic  in  the  x and  y directions.  The  peri- 
odicity in  these  directions  is  characterized  by  the  integers  L and  L , 

x y 

respectively.  Thus,  for  any  function  F which  depends  upon  the  velocities 
and  displacements  of  the  atoms  in  the  lattice  we  have 


; i 


F(x+fcLxao,  y+mLyao,z)  = F(x,y,z) 


(2.1) 


where  l and  m are  arbitrary  integers  and  aQ  is  the  lattice  constant  or 
cube  edge  of  the  conventional  cell. 

The  atoms  will  be  assumed  to  interact  via  a Morse-type  interatomic 
potential  and  therefore  the  Hamiltonian  of  the  lattice  can  be  written 


H.i/2  m .1,2 1 r/*°i?i.AB|-i,.lr  . 

i,o  1,0  i,o  L 


(2.2) 


In  writing  Eq.  (2.2)  we  have  adopted  the  convention  whereby  the  notation 

(i,a)  denotes  the  a1  particle  in  the  i plane  normal  to  the  z axis. 
These  planes  are  numbered  consecutively  beginning  with  the  first  located 
at  z=0.  Any  convenient  labeling  scheme  may  be  used  to  number  atoms 
within  a given  plane,  but  the  particular  convention  used  is  irrelevant 
to  further  discussion  here.  This  labeling  convention  is  convenient  for 
some  of  our  later  discussions. 


In  Eq.  (2.2)  all  quantities  have  been  made  nondimensional . H 
represents,  of  course,  the  total  energy  in  the  lattice  and  has  been 
normalized  by  D,  the  dissociation  energy  of  a single,  isolated  atom  pair; 

v.  is  the  velocity  of  the  (i,a)T  atom,  normalized  by  /D/m,  where  m 

is  the  atomic  mass;  r.  is  the  position  vector  of  the  (i,a)  atom, 

normalized  by  aQ;  Aq  is  the  lattice  constant,  normalized  by  rQ,  the 

separation  of  an  isolated  atom  pair  at  minimum  potential;  and  R is  a 
dimensionless  parameter  representing  the  degree  of  nonlinearity  in  the 
Morse  potential.  The  sum  over  (i,o)  runs  over  all  atoms  in  the  lattice 
and  that  over  (j,8)  is  taken  over  all  atoms  in  the  vicinity  of  the 

(i,a)th  for  which  the  potential  interaction  is  appreciable. 

The  equation  of  motion  satisfied  by  the  (i.a)^  atom  can  be  found 
in  a straightforward  manner  from  Eq.  (2.2)  and  the  result  is 


r.  =2RA  ?. 
i,a  o i,a 


(2.3) 


where  is  the  nondimensional  force  (normalized  by  2 RD/r  ) exerted 

o.,a  th  o' 

on  the  (i,a)  atom  by  the  remaining  atoms  in  the  lattice.  Explicitly, 

is  given  by 


^ » i, 

■ f . 

fit,:'-  >Ar\  V-  S*/ 


F. 

l ,a 


(2.4) 


and  each  dot  represents  differentiation  with  respect  to  the  dimension- 
less time  t,  i.e.,  the  real  time  normalized  by  ►'m/D  aQ. 

In  all  our  calculations  we  will  be  concerned  with  solving  Eqs. 

(2.3)  numerically  to  determine  the  temporal  evolution  of  the  position 
and  velocity  of  each  atom  in  the  lattice  subject  to  some  specific  set 
of  initial  and  boundary  conditions.  From  the  solution  of  these 
equations,  it  is  then  possible  to  infer  all  information  concerning  the 
response  of  the  lattice  to  any  excitation.  The  procedure  for  solving 

Eqs.  (2.3)  is  to  employ  a fourth-order  Runge-Kutta  technique  . The 
details  of  the  method  of  solution,  as  well  as  a listing  of  the  appro- 
priate computer  program,  will  be  presented  elsewhere16  and  will  not  be 
discussed  further  here. 

III.  NUMERICAL  RESULTS  FOR  DISCRETE  LATTICE 

In  all  solutions  of  Eqs.  (2.3)  we  have,  unless  otherwise  stated, 
chosen  the  anharmonicity  factor  R to  be  6.29,  a value  which  leads  to  a 

reasonable  representation  for  a lattice  of  Nickel  atoms1^.  Furthermore, 
we  have  assumed  that  only  atoms  which  were  separated  by  a distance  of 
unity  or  less  (real  distance  normalized  by  aQ)  exerted  an  appreciable 

force  on  one  another.  This  assumption  is  equivalent  to  assuming  that, 
in  the  equilibrium  lattice,  only  an  atom's  first-  and  second-nearest 
neighbors  contribute  significantly  to  its  potential  interaction.  The 
assumption  was  found  to  be  reasonable  for  the  currently  used  value  of 

R.  The  lattice  constant  A.  was  then  calculated16  so  as  to  minimize  the 
potential  and  found  to  be  1.4034. 


15.  B.  Carnahan,  H.A.  Luther,  and  J.O.  Wilkes,  Applied  Numerical 
Methods  (Wiley,  New  York,  1969),  Chap.  6. 

16.  j.d.  Powell  and  J.H.  Batteh,  "Shock  Propagation  in  the  Three-Dimen- 
sional Lattice.  II.  Method  of  Calculation,"  (to  be  published). 

17.  F.  Milstein,  "Applicability  of  Exponentially  Attractive  and  Repul- 
sive Interatomic  Potential  Functions  in  the  Description  of  Cubic 
Crystals,"  J.  Appl.  Phys.  44,  3825  (1973). 
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A. 


Generation  and  Collision  of  Solitary  Waves. 


In  previous  work  we  have  used  the  model  described  in  this  paper  to 
study  shock  propagation  in  solids.  In  those  studies,  we  observed  that, 
when  the  atoms  of  the  lattice  are  initially  at  rest  in  their  equilibrium 
positions,  compression  along  a crystalline  axis  produces  a sequence  of 
solitary  waves  which  propagate  into  the  lattice.  The  solitary-wave 
profiles  used  in  the  computer  experiments  which  are  described  in  this 
paper  were  obtained  by  isolating  a single  soli'  y wave  from  this 
sequence. 


An  example  of  the  shock-wave  calculation  is  shown  in  Figure  1. 
For  this  case,  the  end-most  plane  of  the  lattice,  located  at  z=0,  was 


driven  at  a nondimensional  velocity  Up=1.0  in  the  z direction.  The 


equations  of  motion  of  the  atoms  in  the  lattice  were  solved  and  the 
velocity-time  trajectories  of  various  planes  parallel  to  the  plane  at 
z=0  were  plotted.  For  this  calculation  and  for  all  others  where  the 
motion  of  each  atom  in  a plane  is  identical,  the  same  results  are 


obtained  regardless  of  L and  L , the  periodicity  of  the  lattice  in 

x y 


the  x and  y directions.  Therefore,  and  L^.  can  both  be  set  to  unity 


so  that  we  need  to  solve  the  equations  of  motion  for  only  two  atoms  in 
each  plane,  one  located  at  the  comer  of  the  plane  and  one  located  at 
the  center. 


In  Figure  1 is  plotted  the  z component  of  the  velocity  as  a function 


of  time  x for  the  40l  and  80L  planes  in  the  lattice.  (The  single  sub- 


script is  used  hereafter  to  refer  to  the  plane  as  a whole.)  In  order 
to  facilitate  comparisons  of  the  graphs,  we  have  plotted  the  velocity 


in  each  case  as  a function  of  t-tq  where  tq  is  the  time  at  which  the 
propagating  disturbance  first  excites  the  plane  in  question.  It  is 


evident  from  the  figure  that  a spectrum  of  solitary  waves  is  evolving 
near  the  front  of  the  disturbance,  just  as  occurs  in  one  dimension. 
Asymptotically,  the  pulses  will  completely  separate,  approach  the  same 
constant  shape,  and  propagate  at  a steady  speed  through  the  lattice. 
Consequently,  a single  solitary  wave  can  be  launched  into  a lattice  by 
driving  the  end-most  plane  with  the  solitary-wave  profile  obtained  from 
the  shock-wave  calculations.  This  is  not  the  only  conceivable  method 
of  generating  solitary  waves  in  the  lattice,  but  is  a reasonable  one. 


In  previous  studies  we  demonstrated  that  solitary  waves  propagat- 
ing in  a one-dimensional,  Morse-potential  lattice  are  stable  to  mutual 
collisions.  That  is,  to  within  the  accuracy  of  our  numerical  data,  two 
solitary  waves  emerged  from  a collision  with  the  same  profile  as  prior 
to  the  collision.  It  is  of  interest  to  determine  whether  similar  effects 
occur  in  the  three-dimensional,  FCC  lattice. 
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We  have  launched  two  solitary  waves  having  equal  but  oppositely 
directed  velocities  from  opposite  erds  of  a lattice  and  allowed  them  to 
approach  one  another  and  eventually  collide.  By  plotting  the  velocity- 
time trajectories  of  various  planes  encountered  by  the  solitary  waves, 
the  effect  of  the  collision  can  be  ascertained.  The  results  of  such  a 
calculation  are  shown  in  Figure  2 for  two  solitary  waves  launched  from 
the  ends  of  a lattice  which  was  forty-eight  planes  long.  The  initial 
solitary-wave  profiles  were  obtained  from  a shock-wave  calculation  with 


U =3.0. 
P 


The  top  figure  shows  the  velocity-time  trajectory  of  the  13 


th 


plane  in  the  lattice  beginning  at  the  time  when  the  solitary  wave  propa- 
gating in  the  positive  z direction  first  encounters  the  plane.  For  the 
time  shown  in  this  figure  the  two  waves  have  not  yet  collided  so  that 
this  profile  corresponds  to  the  initial  solitary-wave  profile.  The 
solitary  wave  depicted  in  the  figure  represents  a rather  strong  distur- 
bance. In  fact,  in  the  neighborhood  of  the  solitary  wave  the  density  in 
the  lattice  is  increased  by  about  40%.  The  two  solitary  waves  encounter 


one  another  in  the  vicinity  of  the  24  plane  in  the  lattice  and  the 
trajectory  of  that  plane  is  shown  in  the  center  of  the  figure  during  the 
time  of  collision.  Finally,  at  a much  later  time,  the  collision  has 
been  completed  and  the  negative- velocity  solitary  wave  has  reached  the 

13th  plane.  Its  trajectory  is  shown  in  the  bottom  of  the  figure. 


It  is  apparent  from  the  figure  that  the  two  solitary  waves  maintain 
their  shapes  only  approximately  after  the  collision.  Evidently,  some 
of  the  energy  which  initially  resided  in  the  two  solitary  waves  now 
exists  in  the  form  of  oscillations  which  are  left  behind  by  the  waves. 
Although  it  is  not  apparent  from  the  figure,  the  amplitude  of  the 
solitary  waves  subsequent  to  the  collision  is  decreased  by  about  3% 
from  the  amplitude  prior  to  the  collision.  We  conclude,  then,  that 
solitary  waves  in  an  FCC  lattice  with  a Morse-potential  interaction 
are  not  stable  to  mutual  collisions. 


In  addition  to  the  results  discussed  above,  we  have  also  observed 
collisions  between  solitary  waves  having  smaller  initial  amplitudes. 

As  the  amplitude  decreases,  the  waves  appear  to  become  more  stable. 

In  fact,  solitary  waves  generated  by  steady  compression  with  1/^=1. 0, 

having  initial  amplitudes  of  about  1.8  (as  compared  with  about  5.4  for 
the  Up=3.0  case),  were  found  to  be  stable  to  within  our  numerical 

accuracy. 

It  might  appear  surprising  that,  although  the  only  motion  in  the 
FCC  lattice  is  one-dimensional,  the  solitary  waves  do  not  appear  to  be 
so  stable  as  in  our  previous  one-dimensional -chain  calculations.  The 
reason  for  the  apparent  anomaly  is  not  completely  clear  but  two  possible 
explanations  may  be  offered.  First,  it  must  be  understood  that,  even 
though  the  motion  in  the  Morse-potential  FCC  lattice  is  planar  and  one 
dimensional,  the  model  is  nevertheless  different  from  a one-dimensional 


12 


— > _ 


chain  with  a Morse-potential  interaction.  Consider,  for  instance,  the 
one-dimensional  chain  with  the  atoms  initially  in  their  equilibrium 
positions.  As  two  neighbors  approach  one  another,  their  force  of 
interaction  increases  monotonically,  reaching  a maximum  as  the  separation 
distance  decreases  to  zero.  In  the  FCC  lattice,  however,  as  two  neigh- 
boring planes  approach  one  another,  the  force  exerted  by  an  atom  on  its 
neighbor  in  an  adjacent  plane  is  not  in  the  same  direction  as  that  of 
the  planar  motion.  Therefore,  as  the  planes  approach  one  another,  their 
force  of  interaction  first  increases,  reaches  a maximum  value,  and  then 
decreases  to  zero  as  the  planes  become  coincident.  (Of  course,  it  would 
be  possible  to  reproduce  the  results  of  the  FCC  calculation  with  a one- 
dimensional -chain  model  by  replacing  the  Morse-potential  interaction 
with  another  potential  defined  so  as  to  give  the  same  force  as  a function 
of  separation  distance  in  the  two  cases.)  Second,  it  is  likely  that  the 
solitary  waves  observed  in  our  previous  one-dimensional  calculation  were 
not  rigorously  solitons  but  only  appeared  to  be  so  to  within  our  numerical 
accuracy.  Thus,  after  many  collisions  we  would  expect  that  some  change 
in  the  solitary-wave  profile  would  be  observed.  In  fact,  we  speculated 
then  that  such  an  effect  might  be  responsible  for  the  apparent  tendency 
for  equilibration  that  was  observed  far  behind  the  shock  front  in  our 
one-dimensional  lattice. 


B.  Stability  of  Solitary  Waves  to  Planar  and  Thermal  Oscillations. 

In  addition  to  investigating  the  effects  of  mutual  collisions  upon 
the  solitary-wave  profiles,  we  have  also  examined  the  effects  of 
relatively  small  planar  oscillations  both  along  and  transverse  to  the 
propagation,  direction.  To  perform  the  calculations  we  launched  a 
solitary  wave  at  the  end  of  the  lattice  and  allowed  it  to  propagate  some 
distance  into  the  interior.  At  some  point  ahead  of  the  solitary  wave, 
a few  planes  in  the  quiescent  lattice  were  displaced  slightly,  released, 
and  allowed  to  oscillate.  The  solitary  wave  then  eventually  passed 
through  the  region  of  the  oscillating  planes  and  emerged  on  the  other  side. 
The  intention  was  to  determine  the  effect  of  the  region  of  planar 
oscillations  on  the  solitary-wave  profile. 

An  example  of  the  effect  of  longitudinal  planar  oscillations  is 
shown  in  Figure  3.  A solitary  wave  having  an  amplitude  of  3.6  (generated 
from  compression  at  11^=2. 0)  was  launched  into  the  lattice.  Just  prior 

to  the  time  the  solitary  wave  reached  the  14^  plane  in  the  lattice, 
planes  14-18  were  uniformly  displaced  a distance  of  0.1  in  the  positive 
z direction  and  allowed  to  oscillate.  In  Figure  3a  is  plotted  the  velocity- 
time trajectory  of  the  6*^  plane  in  the  lattice  which  shows  the  unper- 
turbed solitary  wave.  Figure  3b,  on  the  other  hand,  shows  the  same 

disturbance  as  it  propagates  past  the  IS1"*1  plane  in  the  lattice  which 
clearly  lies  within  the  region  of  longitudinal  oscillations.  The  shape 
of  the  original  solitary  wave  is  obviously  distorted  as  it  propagates 


I 


14 


Figure  3.  Effects  of  longitudinal  planar  oscillations  on  solitary-wave 
profile,  (a),  (b) , and  (c)  represent  solitary  wave  before, 
during,  and  after  traversing  the  oscillatory  region, 
respectively. 


through  the  region.  Finally,  in  Fig.  3c  we  have  plotted  the  trajectory  • 

of  the  67*^  plane  in  the  lattice  at  the  time  when  the  original  distur- 
bance first  reaches  it.  By  this  time  the  original  solitary  wave  has 
completely  traversed  the  region  of  longitudinal  oscillations.  Though 
not  shown,  the  trajectories  of  some  later  planes  have  also  been  plotted 
to  demonstrate  that  the  shape  of  the  emerging  solitary  wave  did  not 
change. 

Comparison  of  Figures  3a  and  3c  indicates  that  the  longitudinal 
planar  oscillations  have  had  an  insignificant  effect  upon  the  original 
solitary-wave  profile.  In  fact,  to  within  our  numerical  accuracy,  the 
profile  was  found  not  to  have  changed  at  all.  These  results  might 
appear  surprising  since  the  pulses  were  found  to  be  unstable  to  mutual 
collisions  which  involve  only  longitudinal  planar  oscillations.  Appar- 
’ ently,  however,  the  oscillations  discussed  above  were  too  small  a 

perturbation  to  make  any  change  in  the  wave  profiles  observable  numeri- 
cally. Perhaps  if  the  oscillatory  region  were  much  longer,  or  the 
oscillations  of  larger  amplitude,  a noticeable  effect  would  have  been 
seen. 

Calculations  identical  to  those  above,  except  that  the  five  planes 
were  displaced  in  the  transverse  (y)  direction,  have  also  been  carried 
out.  It  is  interesting  to  note  that  a displacement  of  the  planes  in 
1 the  transverse  direction  will  give  rise  to  both  longitudinal  and  transverse 

oscillations,  whereas  displacement  in  the  longitudinal  direction  produces 
only  longitudinal  oscillations.  (This  point  is  discussed  further  later.) 

Consequently,  we  cannot  examine  the  stability  of  the  solitary  waves  to 
completely  transverse  oscillations,  but  only  to  oscillations  which  contain 
a mixture  of  both.  Furthermore,  although  the  planes  were  displaced  by  the 
same  amount  in  the  two  calculations,  the  change  in  the  energy  in  the 
lattice  due  to  the  transverse  displacement  is  approximately  a factor  of 
two  less  than  that  due  to  the  longitudinal  displacement. 

] 

The  results  of  propagating  the  solitary  wave  through  the  region 
j < containing  transverse  planar  oscillations  are  shown  in  Figure  4.  Again, 

in  Figure  4a  is  shown  the  unperturbed  solitary  wave  and  in  Figure  4b 
the  disturbance  during  the  time  it  traverses  the  region  of  planar  oscil- 
lations. Figure  4c  represents  the  resulting  disturbance  which  emerges 
from  the  oscillatory  region.  As  can  be  seen  from  the  figure,  the  trans- 
verse oscillations  have  significantly  affected  the  size  of  the  original 
l solitary  wave.  In  fact,  the  amplitude  of  the  pulse,  which  was  about 

3.60  prior  to  traversing  the  oscillatory  region  has  been  reduced  to 
about  3.12.  Again,  in  order  to  be  certain  that  the  emerging  pulse  was 
indeed  a solitary  wave,  we  have  followed  its  propagation  farther  than 

the  70**1  plane  into  the  lattice  and  observed  no  change  in  shape. 


Finally,  we  have  allowed  the  solitary  wave  to  propagate  through  a 
region  which  contained  random,  thermal  oscillations.  The  length  of  the 


region  was  the  same  as  for  the  previous  calculations  and  the  cross  section 

contained  32  atoms  (L  =L  =4) . The  thermal  energy  per  particle  was  the 

x y 

same  as  the  energy  per  particle  associated  with  the  mixture  of  trans- 
verse and  longitudinal  planar  oscillations.  Again  some  decay  of  the 
initial  solitary-wave  amplitude  was  observed.  We  should  point  out, 
however,  that  the  finite  size  of  the  lattice  in  the  transverse  directions 
unavoidably  gives  rise  to  some  planar  oscillations.  These  planar  oscil- 
lations, which  are  unexpected  in  macroscopic,  equilibrated  crystals,  no 
doubt  accentuate  the  decay  of  the  solitary-wave  profiles.  Nevertheless, 
we  expect  some  decay  of  the  profiles  even  in  the  absence  of  planar 
oscillations  although  the  decay  should  be  slower  than  that  observed  here. 
Unfortunately,  capability  of  treating  only  small  systems  is  a fundamental 
limitation  of  computer  molecular  dynamics. 

C.  Coupled  Solitary  Waves. 

The  instability  of  the  original  solitary  wave  to  transverse  planar 
oscillations  can  be  explained  in  part  by  the  generation  of  coupled 
longitudinal  and  transverse  solitary  waves.  We  have  observed  these 
waves  in  the  numerical  data  from  which  Figure  4 was  plotted.  Specifi- 
cally, it  was  found  that  there  existed  a solitary  wave  which  propagated 
in  phase  with  the  emerging  longitudinal  wave  but  which  produced  a dis- 
turbance in  the  transverse  (y)  direction. 

The  effect  is  demonstrated  in  Figure  5 in  which  we  have  plotted  the 

velocity-time  trajectory  of  the  70^  plane.  The  upper  part  of  the  graph 
is  identical  to  Figure  4c  and  represents  the  emerging  longitudinal  pulse. 
On  the  lower  part  of  the  graph  is  plotted  the  y component  of  the  velocity 
of  the  same  plane  beginning  at  the  same  time.  As  can  be  seen  from  the 
figure,  the  transverse  solitary  wave  has  somewhat  smaller  amplitude 
than  the  longitudinal  wave  and  propagates  in  phase  with  it.  We  have 
also  computed  the  energy  of  the  initial-wave  profile  and  compared  it 
with  that  of  the  final  coupled  waves.  The  final  energy  was  found  to  be 
about  10%  less  in  the  coupled  configuration,  suggesting  that  part  of 
the  energy  in  the  initial  wave  is  given  up  to  thermal  oscillations. 

We  have  also  observed  coupled  solitary  waves  in  our  study  of  the 
effects  of  thermal  oscillations  upon  the  wave  profile.  In  that  case 
transverse  solitary  waves  were  found  to  propagate  in  both  transverse 
directions.  In  all  likelihood  they  arose  from  the  planar  oscillations 
generated  by  the  finite  size  of  the  crystal  discussed  earlier.  In  the 
following  section  we  will  see  that  the  continuum  limit  of  the  equations 
of  motion  predicts  the  existence  of  these  coupled  solutions  and  we  will 
discuss  their  properties  in  greater  detail. 


IV.  CONTINUUM  EQUATIONS 
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In  this  section  we  will  derive,  interpret,  and  present  some  special 
solutions  to  the  continuum  equations  for  planar  oscillations  in  the 
FCC  lattice.  This  limit  can  be  expected  to  be  valid  whenever  the 
excitations  which  the  equations  describe  have  wavelengths  much  lbnger 
than  the  interatomic  spacing.  The  purpose  of  obtaining  and  solving  the 
equations  is  to  compare  the  results  with  those  obtained  previously  for 
the  discrete  lattice.  As  will  be  seen,  the  relatively  simple  continuum 
equations  predict  results  at  least  qualitatively  similar  to  those 
obtained  in  the  last  section. 


A.  Derivation  of  Continuum  Equations  for  Planar  Oscillations. 

In  order  to  make  the  calculations  as  simple  as  possible  we  will 
assume  that  only  nearest-neighbor  interactions  are  significant.  Includ- 
ing more  distant  neighbors  is  not  expected  to  affect  qualitatively  the 
nature  of  the  results.  This  is  especially  true  for  our  case  where  the 
anharmonicity  factor  is  rather  large.  Furthermore,  since  we  are 
concerned  only  with  planar  oscillations,  each  atom  in  a plane  normal 
to  the  z direction  will  be  assumed  to  have  the  same  velocity  and 
displacement.  The  velocities  and  displacements  may  have  a y component 
(transverse)  and  a z component  (longitudinal)  but,  for  simplicity,  the 
x components  have  been  set  equal  to  zero  throughout. 


We  are  interested  in  solving  Eq.  (2.3),  viz., 
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in  the  limit  in  which  the  displacements  of  all  particles  from  their 
equilibrium  positions  are  small.  The  sum  over  j and  6 now  runs  over  the 
twelve  nearest  neighbors  to  particle  (k,a)  which  are  shown  in  Figure  6. 


For  the  case  of  planar  oscillations,  the  displacement  of  particle 
(k,a)  from  its  equilibrium  position  is  identical  to  the  displacement  of 
the  entire  plane  k from  its  equilibrium  position.  Therefore,  we  can 
write 
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where  the  superscript  o denotes  the  position  vector  to  the  equilibrium 


position  of  the  particle,  and  the  vectors  §.  and  5.  are  the  displacements 

J K 


of  planes  j and  k,  respectively,  from  their  equilibrium  positions.  The 
position  vectors  joining  the  equilibrium  site  of  atom  (k,a)  with  those 
of  its  twelve  nearest  neighbors  are  given  in  Table  I. 


If  we  now  substitute  Eq.  (4.2)  into  Eq.  (4.1),  expand  the  resulting 
equation,  and  retain  only  terms  through  second  order  in  S,  we  obtain 
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TABLE  I.  Position  Vectors  from  Lattice  Site  (k,a)  to  Neighboring 
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Lattice  Sites.  i,j,  and  k are  unit  vectors  in  the  three 
Cartesian  directions. 
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In  deriving  Eq.  (4.3)  we  have  noted  that,  for  nearest-neighbor 
interactions  only,  the  lattice  constant  Re  is  given  by  /2. 
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We  now  perform  the  sum  in  Eq.  (4.3)  using  the  values  of  the  equi- 
librium position  vectors  given  in  Table  I.  After  much  tedious  algebra, 
we  obtain  the  following  equations  for  the  transverse  and  longitudinal 

th 


displacements  of  the  k plane: 
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In  Eqs.  (4.4),  the  subscripts  y and  z denote  the  components  of  the 
displacement  in  the  y and  z directions,  respectively.  Since  we  are 
interested  in  the  case  in  which  the  wavelength  of  the  disturbance  is 

much  larger  than  the  lattice  spacing,  we  can  expand  and  5^  in 

a Taylor  series  about  ^ with  the  result  that 
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Substitution  of  Eq.  (4.5)  into  Eqs.  (4.4)  then  produces 
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where  we  have  now  dropped  the  subscript  k.  Terms  through  fourth  order 
have  been  retained  in  obtaining  Eq.  (4.6). 
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Equations  (4.6)  describe  the  propagation  of  planar  disturbances 
in  the  FCC  lattice  in  the  continuum  limit.  Although  the  result  has  been 
obtained  for  the  Morse  potential,  a similar  result  can  obviously  be 
obtained  for  any  interatomic  potential  whose  associated  forces  are 
expanded  to  second  order  as  in  Eq.  (4.3).  If  we  neglect  third-  and 
fourth-order  terms  (retain  only  the  first  term  on  the  right-hand  side 
of  the  equations)  we  obtain  simply  the  linear  wave  equations.  Thus,  the 
longitudinal  and  transverse  sound  speeds  are  given  in  our  normalization 
by 


and 


= 2 /2  R 
C = 2R  . 


(4.7) 


That  the  ratio  of  these  two  velocities  is  given  by  /2  can  be  inferred 
from  elastic-constant  data  and  the  well-known  fact  that  for  cubic 

crystals*® 


c_i  (c_n\ 
Ct  \ ^44  / 


1/2 


The  ratio  of  the  elastic  constants  C 
two  in  the  limit  of  nearest-neighbor  interactions  only. 


has  been  shown 
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(4.8) 


to  approach 


The  more  interesting  effects  in  Eqs.  (4.6),  however,  are  contained 
in  the  higher-order  terms.  The  third-order  terms  represent,  to  lowest 
order,  the  nonlinear  effects  of  the  potential  whereas  the  fourth-order 
terms  (linear)  account  for  the  dispersive  nature  of  the  lattice.  It  is 
also  clear  that  if  these  higher-order  terms  are  retained,  Eq.  (4.6)  pre- 
dicts that  a transverse  disturbance  cannot  propagate  in  the  absence  of 


a longitudinal  disturbance 
only  trivial,  nonoscillatory  solutions  for  S 


Thus,  setting  Sz=0  in  Eq.  (4.6b)  leads  to 


If,  however,  one  sets 


S =0,  a solution  for  S can  be  found  from  Eq.  (4.6b).  Our  discrete- 

y z 

lattice  results  have  verified  this  effect  as  we  pointed  out  in  Sec.  III. 
B.  In  that  discussion  we  noted  that  displacing  planes  of  atoms  in  the 
transverse  direction  produced  both  longitudinal  and  transverse  oscilla- 
tions, but  a displacement  in  the  longitudinal  direction  produced  only 
longitudinal  oscillations. 


18.  W.P.  Mason,  "Acoustic  Properties  of  Solids,"  in  American  Institute 
of  Physics  Handbook,  edited  by  D.E.  Gray  (McGraw-Hill,  New  York, 
1957),  p.  3-82. 
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Equations  (4.6)  are  obviously  difficult  to  solve  and  we  will  be 
interested  only  in  obtaining  approximate  special  solutions.  Our  primary 
interest  will  be  in  using  the  equations  to  predict  analytically  the 
longitudinal  solitary-wave  profile  and  in  further  studying  the  properties 
of  the  coupled  solitary  waves  discovered  in  Sec.  III. 

The  equations  can  be  somewhat  simplified  by  noting  that  solitary 
waves  represent  steady,  travelling-wave  solutions  to  the  equations  of 
motion.  Thus,  we  assume  solutions  of  the  form 


sy  = sy(k-CT)  = sy(e) 
sz  = sz(k-CT)  = sz(o 


(4.9) 


where  C is  the  propagation  velocity  of  the  resulting  wave  form,  and 
substitute  into  Eqs.  (4.6).  If  we  define  the  components  of  the  planar 
velocity  as 

3S 

(4.10) 


Eqs.  (4.6)  become 


u"  = au  - 4Buv 
v* ' = yv  - 5v2-8u2 


(4.11a) 

(4.11b) 


The  primes  denote  differentiation  with  respect  to  £ and  the  constants 
are  defined  as  follows: 

a = 12(C2/C2-1) 


6 = 3 (3R- 1) /C 

y = 12 (C2/C2-l) 


(4.12) 


6 = 18  (R-l)/C  . 

In  the  following  two  sections  we  obtain  and  discuss  some  solutions  to 
Eqs.  (4.11). 


taw#; 


B.  Longitudinal  Solitary  Waves. 


The  simplest  solution  of  Eqs.  (4.11)  results  whenever  no  transverse 
disturbance  exists  so  that  u can  be  set  equal  to  zero  throughout.  We 
then  have 

v//  =*  yv  - 6v2  . (4.13) 

If  Eq.  (4.13)  is  multiplied  by  V , it  can  be  integrated  twice  provided 
that  we  require  that  the  function  and  its  derivatives  vanish  at  infinity. 
Performing  the  calculation  we  find 

v = sech2  (^  5)  • (4.14) 

In  obtaining  Eq.  (4.14)  <s  have  assumed  y>0;  the  solution  for  y<0  is 
oscillatory,  does  not  vanish  at  infinity,  and  will  not  be  considered 
here.  Analytic  approximations  for  solitary-wave  profiles  similar  to 

19 

this  result  have  been  discussed  in  the  literature  . 

In  order  to  compare  the  wave  profile  predicted  by  Eq.  (4.14)  with 
our  numerical  data  for  the  discrete  lattice,  we  plotted  the  velocity- 
time trajectory  of  a plane  of  atoms  in  the  lattice  as  a solitary  wave 
traversed  it.  The  propagation  velocity,  C,  of  the  solitary  wave  was 
obtained  from  the  computer  data  and  substituted  into  Eqs.  (4.12).  The 
resulting  constants  were  then  used  to  calculate  the  profile  in  Eq.  (4.14). 
Results  of  the  calculation  are  shown  in  Fig.  7 in  which  are  plotted  the 
numerical  and  analytic  wave  profiles  as  a function  of  the  parameter  £/C. 
The  two  profiles  are  in  reasonably  good  agreement.  We  have  also  compared 
a number  of  other  profiles  for  both  higher-  and  lower-  amplitude  solitary 
waves.  However,  as  the  amplitude  of  the  wave  increases,  its  width 
decreases.  The  discrepancy  between  analytic  and  numerical  results  then 
increases  owing  to  the  inadequacy  of  the  continuum  approximation. 

C.  Coupled  Solitary  Waves. 

We  now  wish  to  determine  whether  we  can  predict  from  Eqs.  (4.11)  the 
coupled  longitudinal  and  transverse  solitary  waves  observed  previously. 

We  begin  by  obtaining  an  approximate  analytic  solution  to  the  equations 
valid  in  the  limit  u -*■  0.  The  technique  may  be  viewed  as  the  first  step 
in  an  iterative  procedure;  in  principle,  it  may  be  repeated  as  often  as 
desirable. 


19.  N.J.  Zabusky,  "A  Synergetic  Approach  to  Problems  of  Nonlinear 
Dispersive  Wave  Propagation  and  Interaction,"  in  Nonlinear 
Partial  Differential  Equations,  edited  by  W.F.  Ames  (Academic, 
New  York,  1967),  p.  223. 
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Figure  7.  Comparison  of  numerical  and  analytic  longitudinal  solitary- 
wave  profiles.  The  dashed  curve  represents  the  numerical 
results  for  the  discrete  lattice;  the  solid  curve  represents 
a plot  of  Eq.  (4.14) . 


Equation  (4.14)  represents  the  solution  of  bq.  (4.11b)  for  u=0. 
Substituting  the  result  into  Eq.  (4.11a)  we  obtain  an  equation  for  the 
first  approximation  to  u,  namely, 


u' ' = au  - sech2  ^ E u . 


(4.15) 


Equation  (4.15)  is  linear  and  second  order  and  is  similar  to  the  time- 

independent  Schroedinger  equation.  Since  a>0,  one  expects20  that 
bounded  solutions  exist  only  for  certain  discrete  values  of  C and  that 
the  corresponding  functions  u vanish  at  infinity. 


To  solve  the  equation  we  make  the  change  of  variable 


/7 

y = sech  £ 


and  substitute  into  Eq.  (4.15)  to  produce 
2 

y 2,,  2.  d u y ,,  _ 2.  du 

4 y (i-y  ) — ♦ j yd-2y  ) 37  = «u  - 
dy  7 


66y  2 

6 y u 


(4.16) 


(4.17) 


Wc  now  assume  that  the  solution  to  Eq.  (4. 17)  can  be  represented  in  the 
form  of  the  series 


v n 

u = l V 

p=o  ^ 


(4.18) 


where  m is  a positive  number,  and  substitute  into  Eq.  (4.17).  Equating 

coefficients  of  y in  the  resulting  expression  then  yields  the  recursion 
relation 


2g.-2)  (m+2ft-l)-; 
y(m+2fc)2-4a 


(4.19) 


By  assumption  uq*0  and  un=0  for  n < 0 so  that  the  denominator  of 
Eq.  (4.19)  must  vanish  for  t=0.  Thus,  we  obtain  the  eigenvalue  equation 


(4.20) 


20.  A.  Messiah,  Quantum  Mechanics  (Wiley,  New  York,  1966),  Vol.  1, 
Chap.  2. 
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which  becomes 


I 

f 


f 

t 

I 


m2-4 

m2-8 


(4.21) 


after  using  Eqs.  (4.12).  This  result  indeed  predicts  that  solutions  to 
Eq.  (4.17)  exist  only  for  certain  well-defined  values  of  C.  Furthermore, 
if  the  series  represented  by  Eq.  (4.18)  is  to  terminate,  the  numerator 
of  Eq.  (4.19)  must  vanish  for  some  value  of  i.  This  condition  implies 


m - ~ (4fc-3)  ♦ [1»16(3R-1)/ (R-l) 1 1/2 
in  — _ 


(4.22) 


where  all  constants  have  been  expressed  in  terms  of  R.  Of  course,  m 
must  be  greater  than  zero  in  order  for  the  solution  to  remain  bounded  at 
infinity. 

As  a specific  example,  let  us  now  evaluate  the  solution  for  the 
case  R=6.29.  In  that  case,  Eq.  (4.22)  predicts  that  the  only  solution 
which  leads  to  m>0  and  y>0  occurs  for  i.=  l and  m=  3 . 2 1 . Substitution  of 
m into  Eq.  (4.21)  then  yields  0=1.660^  and,  from  this  value  of  C,  the 

remaining  constants  can  be  calculated  from  Eq.  (4.12).  Only  the  zeroth 
term  then  survives  in  the  expansion  of  Eq.  (4.18)  and  we  have 

u = uq  sech3,21(2.3£) . (4.23) 

The  corresponding  expression  for  v,  obtained  from  Eq.  (4.14),  is 

v = 9.8  sech2  (2.3£).  (4.24) 

In  principle,  Eq.  (4.23)  could  now  be  used  in  Eq.  (4.11b)  to  obtain  a 
second  approximation  to  v. 

It  is  important  to  emphasize  that  Eqs.  (4.23)  and  (4.24)  represent 
only  a rather  crude  approximation  to  Eq.  (4.11)  which  can  be  expected 
to  be  valid  only  as  u tends  to  zero.  The  solution  predicts,  for  instance, 
that  the  amplitude  of  the  longitudinal  solitary  wave  in  the  coupled-wave 
profile  is  identical  to  that  for  an  isolated  longitudinal  wave  having 
the  same  value  of  C;  actually  the  amplitude  is  diminished  from  that  of 
the  isolated  wave.  Furthermore,  the  value  of  u cannot  be  obtained  in 

the  lowest-order  solutions  to  the  equations  and  this  parameter  must  be 
fit  to  the  data.  Finally,  the  lowest-order  approximation  yields  only 
one  acceptable  eigenvalue  and,  thus,  one  solution  to  the  equations. 

Other  solutions  are,  however,  possible  as  u increases  from  zero  and 
these  are  no  doubt  predicted  by  higher-order  approximations.  Despite 
these  shortcomings,  the  analytic  solution  is  nonetheless  valuable 
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because  it  does  predict  that  coupled-wave  solutions  do  exist  and,  as 
we  shall  see,  predicts  the  shapes  of  the  wave  profiles  rather  well  for 


Because  of  the  limitations  of  the  analytic  solution  we  have  also 

obtained  numerical  solutions  of  Eqs.  (4.11).  Our  efforts  have  been 

confined  to  obtaining  solutions  which  propagate  in  phase,  that  is,  to 

solutions  which  reach  their  maximum  at  the  same  value  of  £.  Only 

solutions  having  this  property  have  been  observed  in  our  studies  of  the 

discrete-lattice  equations.  The  numerical  procedure  was  as  follows: 

A value  for  v was  assumed  (arbitrarily  chosen  at  £=0)  and  u w 
max  v ' max 

obtained  from  the  differential  equations.  Specifically,  if  we  multiply 

Eq.  (4.11a)  by  u ' , Eq.  (4.11b)  by  2v'  , and  add,  the  resulting  expression 

is  in  the  form  of  an  exact  differential.  Integrating  and  evaluating 

the  resulting  equation  at  £=0  yields 


2 6v  /3-y 

max 

a/2-2Bv 


(4.25) 


These  values  of  u and  v were  used  as  initial  conditions  and  a 
max  max 

fourth-erder  Runge-Kutta  scheme  employed  to  solve  Eqs.  (4.11).  As 
expected  from  the  analytic  result,  the  numerical  solutions  were  found 
to  diverge  as  Id-*00  unless  the  appropriate  value  of  C.  found  by  trial 
and  error,  was  used. 

We  have  performed  a number  of  numerical  solutions  of  Eqs.  (4.11) 

in  the  manner  discussed  assuming  various  values  of  v .In  no  case 

max 

were  we  able  to  obtain  convergent  solutions  for  values  of  v less 

& max 

than  or  equal  to  the  value  predicted  by  the  analytic  technique  in  lowest 

order,  namely,  v =9.8.  The  result  suggests  that  there  exists  a threshold, 
1 max 

occurring  at  the  eigenvalue  C=1.66Cjl,  below  which  the  coupled  solitary 

waves  cannot  propagate  in  phase.  Unfortunately,  the  threshold  is 
sufficiently  high  that  one  cannot  expect  that  the  continuum  equations  are 
a reasonable  quantitative  approximation  to  the  discrete-lattice  equations 
and  a direct  comparison  of  the  results  is  not  possible.  Nonetheless,  it 
is  interesting  that  these  equations  predict  both  the  existence  of  the 
coupled  waves  as  well  as  the  fact  that  they  obey  a threshold. 


An  example  of  a numerical  solution  is  shown  by  the  solid  curve  in 
Figure  8.  The  solution  was  obtained  by  assuming  a value  of  vmax  given 

by  10.24  and  was  found  to  converge  provided  C was  given  by  29.9  or 
1.68  C^.  A plot  of  u obtained  from  the  resulting  solution  is  shown  in 

the  lower  half  of  the  figure;  umax  was  found  to  be  approximately  1.5 

as  can  be  verified  from  Eq.  (4.25). 
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Since  u is  considerably 
solution  represented  by  Eqs. 
approximation  to  the  numeric 
case  as  can  be  seen  from  the 
sents  the  solution  given  by 

Eq.  (4.23)  equal  to  1.5  and 
half  of  the  graph  along  with 
were  found  to  be  coincident 
plot  the  data.  Consequently 
a reasonable  approximation  t 


values  of  v 


near  the  thre 


less  than  v one  expects  that  the  analytic 
(4.23)  and  (4.24)  might  well  be  a reasonable 
al  result.  This  is  in  fact  found  to  be  the 
top  graph  in  which  the  dashed  curve  repre- 
liq.  (4.24).  Furthermore,  when  we  set  uQ  in 

attempted  to  plot  the  results  on  the  lower 
the  numerical  solution  for  u,  the  graphs 
to  within  the  accuracy  with  which  we  could 
, we  conclude  that  the  analytic  solution  is 
• the  coupled  solitary-wave  profile  for 
shold  value  of  9.8. 


Equations  (4.11)  merit  further  study.  We  have  confined  our  attention 
in  this  investigation  only  to  solutions  which  have  their  maximum  at  the 
same  value  of  £.  We  have  been  unable  to  prove  from  the  differential 
equations,  however,  that  all  solutions  which  vanish  at  infinity  have 
this  property.  Thus,  although  we  have  not  observed  it  in  solutions  to 
the  discrete-lattice  equations,  there  remains  the  possibility  that  coupled 
waves  exist  which  propagate  out  of  phase.  It  is  interesting  to  note 
that  if  such  solutions  exist  they  too  obey  a threshold  effect,  although 
not  necessarily  the  same  one  discussed  previously.  In  fact,  it  is 
shown  in  the  Appendix  that  for  all  bounded  solutions  to  Eqs.  (4.11), 

C must  obey  the  condition 

/, o \ 1/2 


'(■*») 


(4.26) 


For  R=6.29,  this  yields  C > 1 . 31C^  which  is  substantially  lower  than 
the  value  C > 1.66  C suggested  for  the  in-phase  solutions. 


V.  SUMMARY  AND  CONCLUSIONS 

We  have  undertaken  both  computer-molecular-dynamic  as  well  as  some 
analytic  studies  of  sol itar> -wave  propagation  in  the  three-dimensional, 
FCC,  Morse-potential  lattice.  It  has  been  found  that  the  lattice  is 
capable  of  supporting  the  propagation  of  solitary  waves  as  in  one- 
dimension.  The  basic  conclusions  reached  regarding  the  properties  of 
the  solitary  waves  are  as  follows: 

1.  Longitudinal  solitary  waves  are  not  stable  to  mutual 
collisions,  the  degree  of  stability  decreasing  as  the 
amplitude  of  the  colliding  solitary  waves  increases. 
Nevertheless,  even  in  a mutual  collision,  which  repre- 
sents a rather  strong  longitudinal  disturbance,  the 
solitary  wave  rctans  its  integrity  fairly  well. 


. >v  ..v”  / +■ 


2.  The  solitary  waves  are  extremely  stable  to  small,  longi- 
tudinal, planar  oscillations. 

3.  Longitudinal  solitary  waves  are  not  stable  to  small, 
transverse,  planar  oscillations.  However,  even  in  this 
case,  most  of  the  energy  associated  with  the  initial 
solitary  wave  remains  localized  in  the  form  of  a 
coupled  longitudinal  and  transverse  solitary  wave. 

4.  Solitary  waves  appear  to  be  more  stable  to  random 
thermal  oscillations  than  to  coherent  planar 
oscillations,  although  it  is  difficult  to  determine 
this  conclusively  from  computer  simulations. 

Our  calculations  suggest  that,  for  a variety  of  perturbations, 
the  energy  initially  associated  with  a solitary  wave  tends  to  remain 
localized  within  the  wave,  although  there  may  be  an  exchange  of  energy 
between  longitudinal  and  transverse  oscillations.  Therefore,  despite 
the  absence  of  total  stability  of  the  solitary  waves,  they  may  none- 
theless be  important  in  three-dimensional  energy-transport  problems. 

For  instance,  computer  simulations  have  shown  that  solitary  waves  are 
generated  whenever  a solid  is  subjected  to  shock  compression  and  they 
will,  no  doubt,  be  generated  in  other  nonequilibrium  problems  as  well. 
Their  fair  degree  of  stability,  then,  insures  that  solitary  waves  will 
substantially  affect,  at  least  initially,  both  the  relaxation  time  and 
the  manner  in  which  thermal  equilibrium  is  re-established  after  a 
disturbance.  Furthermore,  for  systems  of  realistic  dimensions  which 
are  initially  in  thermal  equilibrium,  the  solitary  waves  are  likely  to 
be  more  stable  than  the  calculations  reported  here  suggest  because  of 
the  absence  of  coherent,  planar  oscillations  in  the  background. 

Many  problems  of  importance  to  ballistics  research  require  an 
understanding  of  the  manner  in  which  energy  is  transported  in  solids. 

For  example,  models  for  the  initiation  of  explosives  or  for  the  defor- 
mation of  materials  under  impulsive  loads  must  incorporate  some  assump- 
tions regarding  the  nature  of  the  energy  transfer  and  relaxation  pro- 
cesses which  occur  following  a disturbance.  If  these  processes  in 
real,  impure  crystals  are,  in  fact,  dominated  by  solitary-wave  be- 
havior, then  it  will  become  necessary  to  include  the  effects  of  the 
solitary  waves  in  the  models.  For  instance,  solitary  waves  may  sig- 
nificantly affect  the  rate  of  chemical  reactions  behind  a shock  wave 
since  the  velocity  distribution  function  is  highly  non -Maxwellian  in 
a region  where  solitary  waves  are  propagating.  Consequently,  future 
research  should  be  directed  towards  determining  the  extent  to  which 
solitary  waves  persist  in  more  realistic  (e.g.,  impure)  crystals  and 
assessing  qualitatively  their  effect  on  the  macroscopic  properties 
of  the  lattice.  Eventually,  it  would  be  desirable  to  develop  a 
methodology  for  incorporating  detailed,  quantitative  studies  of  solitary- 
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wave  behavior  into  the  macroscopic  description  of  energy  transfer  in 
real  solids.  This  goal,  however,  will  require  a significant  advance 
of  the  state-of-the-art  both  in  computer-molecular-dynamic  techniques 
and  in  our  knowledge  of  atomic  and  molecular  interactions. 
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The  purpose  of  this  appendix  is  to  prove  that  for  solutions  of 
Eqs.  (4.11)  which  vanish  at  infinity,  the  eigenvalues  obey  a threshold 
condition.  That  is,  we  demonstrate  that  for  solutions  of  the  equations 
to  exist  we  must  have 


C / 2 > 2R- 

' Cl  1+R 


Consider  Eqs.  (4.11),  viz., 

u77  = au  - 4 Buv 


(A.l) 


(A. 2) 


v7 7 ~ yv  - 6v2  - Bu 


If  the  solution  for  u is  to  vanish  at  infinity,  then  it  must  have  a 
maximum  in  a region  where  u is  positive  and/or  a minimum  in  a region 
where  u is  negative.  The  final  result  is  independent  of  which  of 
these  situations  exists,  so  let  us  assume  that  u has  a maximum  value 
u*  at  £ = £*,  and  that  u*>0.  Furthermore,  let  us  denote  the  value  of 
v at  £*  by  v*.  If  u is  to  be  a maximum  at  C = £*,  its  second  derivative 
evaluated  at  must  be  negative.  Consequently,  from  Eq.  (A. 2)  we  have 


and,  since  u*  > 0, 


(a-4Bv*)u*  < 0 


4B 


(A. 3) 


(A. 4) 


If  we  multiply  the  first  of  Eqs.  (A. 2)  by  u7  and  the  second  by  2V7  , 
add,  and  integrate,  we  obtain 


1/2 (u7  ) + (v7  )‘  = ♦ Yv‘  - 6vJ  - 2Bu"vX)  . 

We  now  apply  relation  (A. 5)  at  £ = £*  to  produce 
1/2  (a-4Bv*)u*2  + v*2(y-2/3Sv*)  > 0 . 

Equations  (A. 4)  and  (A. 6)  imply  that  v*  must  obey  the  relation 


< v*  < ^ 
46  26 
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(A. 5) 


(A. 6) 


(A.  7) 


1 


In  order  for  a solution  to  exist,  then,  we  must  have  that 


a 

48 


(A.  8) 


Expressing  the  constants  in  terms  of  C and  R [see  Eqs.  (4.12)]  yields 
finally 


(A. 9) 
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